home *** CD-ROM | disk | FTP | other *** search
/ Ham Radio 2000 #2 / Ham Radio 2000 - Volume 2.iso / HAMV2 / ANTENNA / YAGIU112 / ZBR_HACK.C < prev   
Encoding:
C/C++ Source or Header  |  1995-08-07  |  1.3 KB  |  72 lines
  1. #include <math.h>
  2. #include "yagi.h"
  3.  
  4. #define ITMAX 100
  5. #define EPS 3.0e-8
  6.  
  7. double  zbrent(double (*error_3dB) (double x),double x1,double x2, double tol)
  8. {
  9.     int iter;
  10.     double a=x1,b=x2,c,d,e,min1,min2;
  11.     double fa=(*error_3dB)(a),fb=(*error_3dB)(b),fc,p,q,r,s,tol1,xm;
  12.     void nrerror();
  13.     /* printf("x1=%lf fa=%lf x2=%lf fb=%lf\n", x1,fa,x2,fb); */
  14.     if (fb*fa > 0.0) nrerror("Root must be bracketed in ZBRENT");
  15.     fc=fb;
  16.     for (iter=1;iter<=ITMAX;iter++) {
  17.         if (fb*fc > 0.0) {
  18.             c=a;
  19.             fc=fa;
  20.             e=d=b-a;
  21.         }
  22.         if (fabs(fc) < fabs(fb)) {
  23.             a=b;
  24.             b=c;
  25.             c=a;
  26.             fa=fb;
  27.             fb=fc;
  28.             fc=fa;
  29.         }
  30.         tol1=2.0*EPS*fabs(b)+0.5*tol;
  31.         xm=0.5*(c-b);
  32.         if (fabs(xm) <= tol1 || fb == 0.0) return b;
  33.         if (fabs(e) >= tol1 && fabs(fa) > fabs(fb)) {
  34.             s=fb/fa;
  35.             if (a == c) {
  36.                 p=2.0*xm*s;
  37.                 q=1.0-s;
  38.             } else {
  39.                 q=fa/fc;
  40.                 r=fb/fc;
  41.                 p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0));
  42.                 q=(q-1.0)*(r-1.0)*(s-1.0);
  43.             }
  44.             if (p > 0.0)  q = -q;
  45.             p=fabs(p);
  46.             min1=3.0*xm*q-fabs(tol1*q);
  47.             min2=fabs(e*q);
  48.             if (2.0*p < (min1 < min2 ? min1 : min2)) {
  49.                 e=d;
  50.                 d=p/q;
  51.             } else {
  52.                 d=xm;
  53.                 e=d;
  54.             }
  55.         } else {
  56.             d=xm;
  57.             e=d;
  58.         }
  59.         a=b;
  60.         fa=fb;
  61.         if (fabs(d) > tol1)
  62.             b += d;
  63.         else
  64.             b += (xm > 0.0 ? fabs(tol1) : -fabs(tol1));
  65.         fb=(*error_3dB)(b);
  66.     }
  67.     nrerror("Maximum number of iterations exceeded in ZBRENT");
  68. }
  69.  
  70. #undef ITMAX
  71. #undef EPS
  72.